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Abstract 

This paper presents the generahzation of weighted distances to modules and their 
computation through the chamfer algorithm on general point-lattices. A first part 
is devoted to the definitions and properties (distance, metric, norm) of weighted 
distances on modules, with a presentation of weight optimization in the general case, 
to get rotation invariant distances. A general formula for the weighted distance on 
any module is presented. The second part of this paper proves that, for any point- 
lattice, the sequential 2-scan chamfer algorithm produces correct distance maps. 
Finally, the definitions and computation of weighted distances are applied to the 
face-centered cubic (FCC) and body-centered cubic (BCC) grids. 

Key words: Weighted distance. Distance transform, Chamfer algorithm, 
Non-standard grids 



1 Introduction 



Given a binary image consisting of object and background grid points, the 
distance transform assigns to each object grid point its distance value to the 
closest background grid point. The distance transform provides valuable infor- 
mation for shape analysis and is widely used in several applications [1], such 
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as, for example, skeleton extraction [2], template matching [3], shape based 
interpolation [4,5] or image registration [6]. 

Computing the distance from each object grid point to each background grid 
point would lead to a far too high computational cost. Numerous authors have 
thus investigated alternative ways of computing distance maps. To do so, they 
use the spatial consistency of a distance map, which allows propagation of 
local information. In this paper, we focus on one such algorithm, the chamfer 
algorithm, which computes the weighted distance transform (WDT). When 
using the chamfer algorithm, only a small neighborhood of each grid point is 
considered. A weight, a local distance, is assigned to each grid point in the 
neighborhood. By propagating the local distances in the two-scan algorithm, 
the correct distance map is obtained. For example, the well-known two dimen- 
sional city-block (L^) and chessboard {L°°) distances can be obtained in this 
way by using unit weights for the neighbors. 

It has been shown that in two dimensions, the hexagonal grid is in many ways 
to prefer over the usual square grid, [7]. For example, since the hexagonal grid 
constitutes the closest sphere packing [8] , only 87% of the number of samples 
needed in the square grid can be used without information loss, [9] . Extending 
this to three dimensions, the face-centered cubic (FCC) grid is the lattice with 
the highest packing density, resulting in that less samples are needed for the 
FCC grid without loosing information, compared to the cubic grid, [9]. The 
reciprocal grid of the FCC grid is the body-centered cubic (BCC) grid. This 
means that if the FCC grid is used in the spatial/frequency domain, then 
the BCC grid is used to represent the image in the frequency/spatial domain. 
The FCC grid has the highest packing density, so this grid is preferably used 
in frequency domain resulting in that only 71% of the samples are needed 
for an image on the BCC grid compared to the cubic grid. If the FCC grid 
is used in spatial domain instead, 77% of the samples can be used without 
loosing information, [9]. This reasoning requires that it is possible to acquire 
images directly to the grids with non-cubic voxels. In theory, any tomography- 
based technique for acquiring volume images (i.e., that gives a sequence of 2D 
slices) such as Magnetic Resonance Imaging (MRI), Computed Tomography 
(CT), Single Photon Emission Computed Tomography (SPECT), or Positron 
Emission Tomography (PET) can be adjusted to work on the FCC or BCC 
grids. 

For SPECT and CT, using the filtered back-projection method or direct 
Fourier methods, the fact that the Fourier transform for the FCC and BCC 
grids (and for two-dimensional planes that correspond to the 2D slices) exists 
is enough. Images on the FCC and BCC grids acquired using the algebraic re- 
construction technique for SPECT or CT is found in [10]. They use spherically- 
symmetric volume elements (blobs), [11]. With the PET-technique, the origin 
of gamma-rays (produced by annihilation of an emitted positron and an elec- 
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tron) are computed. The technique does not depend on the underlying grid 
except in the digitization step, which just as well can be performed on the 
FCC and BCC grids. 

Dealing with digital images means dealing with discrete grids. Indeed, digi- 
tal images are discrete representation of the continuous world. Two kind of 
approaches can be adopted to deal with discrete grids. The first one applies 
continuous methods to the discrete grid. This is the case of the Euclidean 
distance transform (EDT) which aims at computing exact Euclidean distance 
on discrete images. The second approach, called discrete geometry, consists in 
developing tools directly devoted to digital images and discrete grids. This is 
the case of the WDT which can, for example, be computed using a two-scan 
algorithm, a chamfer algorithm in the usual square grid [12]. 

Since computing the WDT is less memory demanding and faster than com- 
puting the EDT, the WDT is preferable in situations where the data set is 
too big to fit into the central memory (such as in [13]) or the computational 
time should be minimized (recent contributions involve speedup of level-set 
methods, [14], and real-time tele-operated force- feedback, [15]). It should be 
mentioned that there are linear-time (or almost linear) algorithms for comput- 
ing the EDT obtained by sequentially propagating vectors, [16], and dimen- 
sionality reduction, [17,18,19]. These algorithms are, however, not as fast as 
the chamfer algorithm since the EDT requires a larger number of arithmetic 
operations. Also, the vector propagation algorithm is never guaranteed to be 
error-free, [20] . The WDT has other properties which makes it useful in appli- 
cations. It consists of only integers without loosing the metric property which 
is not the case for the EDT (real- valued distances) or the squared EDT (does 
not satisfy the triangular inequality). This leads to nice properties, e.g., when 
extracting the centers of maximal balls for computing skeletons which is fast 
and simple for the most common WDTs, [21,1], in contrast to the case for 
EDT, [22]. Observe that, when using large neighborhoods in the WDT, then 
the situation can be complex also for WDTs, [23]. In general, the weighted 
distances are well-suited for morphological operations, [24]. 

To allow image processing on images on the non-standard grids, algorithms 
must be designed to work directly on these grids. Weighted distances and the 
chamfer algorithm has been applied to non-standard grids such as the two- 
dimensional hexagonal grid, [25,26,27] and the FCC and BCC grids, [27]. In 
this paper, the weighted distance and the chamfer algorithm are considered 
on modules. This gives a very general framework in which Z", the hexagonal, 
FCC, and BCC grids are all considered in parallel. The general framework 
implies that the theory also applies to the higher- dimensional generalizations 
of these grids. 

This paper is organized as follows: we first summarize weighted distance defi- 
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nitions and properties that can be found in the hterature and generahze them 
to sub-modules of M". Then we exhibit conditions for the chamfer algorithm 
to work on our framework. Finally, we propose examples of chamfer masks for 
the body-centered cubic (BCC) grid and the face-centered cubic (FCC) grid. 



2 Definitions, notations and properties 

This section generalizes definitions and weighted distance properties found 
in literature to our general framework of modules. We denote the set of real 
numbers, the set of integers, and the set of natural numbers, with M, Z, and 
N, respectively. 

2.1 General framework: Module 

Let {Q, +) be an Abelian group {Q could be, for example, M" or Z" for n G 
N\{0}). To define weighted distances, we need not only an internal operator 
(here +), but also an external operation: a ■ p with a G M, Z, or N and p G ^. 
Vector spaces (as for example (]R",]R, +, x)) handle external operations, but 
are too restrictive ((Z'^,Z, +, x), for example, is not a vector space). To be 
general enough, we use modules. 

Definition 2.1 (Module) LetTZ be a commutative ring (for example (Z, -|-, x)^ 
with two neutral elements and 1. A set Q is called a module on 7^ (or TZ- 
module) if Q has a commutative group operation +, an external law ■ , and 
satisfies the following properties: 

(identity) Vp G ^ 1 ■ p = p 

(associativity) Vp G Q, Va, P eTZ, a ■ (/5 ■ p) = (a x /3) ■ p 

(scalar distributivity) Vp G Q , Va, P eTZ, (a + /3)-p = a- p + /5-p 

(vectorial distributivity) Vp, q G ^, Va G 7?., a ■ (p + q) = a ■ p + a ■ q 



We now consider sub-rings (7^,+, x) of (M, +, x). Given an Abelian group 
{G, +) we consider the 7?.-module {Q, 71, +, ■). 

Remark 2.1 The main difference between a module and a vector space is 
the non-invertibility (with respect to the external law) of the elements of its 
associated ring TZ (e.g. 2 G Z but 1/2 ^ Z). A basis of a module Q of dimen- 
sion n is a family of n independent vectors (vi)i=i..n (Va G 7^", aj ■ vf = 
-v^ Vi G [l..n], Q!j = 0). But a linearly independent family of n vectors may 
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not be a basis of Q. For example, ^(1,0), (0,2)^ is a family of two indepen- 
dent vectors of but is not a basis of ((1, 1) G can not be reached by 
a linear combination of (1, 0) and (0, 2) with coefficients taken in the ring Z). ♦ 

Definition 2.2 (Distance) A distance on a group Q, having values in TZ, 
called {d, Q, TZ) is a function d:QxQ\-^7l which satisfies the following prop- 
erties: 

(positive) Vp, q G Q d{p,q) > 
(definite) Vp, qE Q d{p, q) = <(=^ p = q 
(symmetric) Vp, qE Q d{p, q) = d{q, p) 

Definition 2.3 (Metric) Given a distance {d,Q,Tl), d is called a metric if 
it also satisfies the following property: 

(triangular inequality)\/p, q,r E Q d{p, q) < d{p, r) + d(r, q) 

Definition 2.4 (Norm) Given a metric {d,Q,TZ), d is called a norm on the 
module {Q, TZ, +, ■) if it also satisfies the following property: 

(positive homogeneity) Vp, q G ^, Va G 7^, d{a ■ p, a ■ q) = |q;| x d{p, q) 

Definition 2.5 (Image) An image is a function I: S ^ TZU oo, where S is 
a finite subset of Q . 

Definition 2.6 (Distance map) Given a binary image I , i.e., /(p) G {0, l},let 
X = {p E S, /(p) = 1} be the foreground and X = {p E S, J(p) = 0} be the 
background. Given a distance {d,Q,7lUoo), the distance map, DMx, of I 
is a grey level image, DMx{p) G TZ, where the value of each point of the 
foreground corresponds to its shortest distance to the background, i.e. 

is — >nuoo 

DMx : <^ _ 

[ p I — > d{p, X) = infqgY d{p, q) 

Given a module (i.e. a set of vectors) [Q, TZ) and an affine space (i.e, a set of 
points in /) Q, we say that Q is equivalent to Q in the the following sense: 

Vp, q e ^, 33? G ^, q = p + 3?; 
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this vector it is denoted = Q ~ P- For every module Q, there exists an 
equivalent affine space Q and vice versa. These spaces always have the same 
dimension. Given a point O of Q, {Q, O) is an affine space with an origin. 
An origin is required to build a basis of Q and allows to define the following 
operations. 

Vp, q e ^, VA e 7^, Ap = O + AOp, p + q = + Op + Oq. 
In the following, we will consider Q and Q as the same set and denote p = Op. 

2.2 Chamfer masks and weighted distances 

We now consider more particularly sub-modules of (M",M), (for example 
(Q", Q) and (Z", Z)) equivalent to the affine space M", with the origin 0(0, 0, 0) 
and the canonical basis {(1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 0, 1)}. For any 
vector or point x G ^, we denote x = {x'^)i=i..n the decomposition of x in this 
canonical basis. 

Among numerous other methods, distance maps can be computed by propa- 
gation of local distances using a chamfer mask. The latter is defined as follows: 

Definition 2.7 (Chamfer mask) A chamfer mask C is a finite set of weighted 
vectors {i^k,Wk)k€[i..m] x 7^} which contains a basis of {Q,TZ) and satis- 
fies the following properties: 

{positive weights) Wk Wk G IR+ and 7^ 
{symmetry) (V, w) E C =^ {—^: w) E C 

Intuitively, we speak about a distance between two points as the length of 
the shortest path between these points. In the case of weighted distances, we 
restrict the possible paths to those allowed by the chamfer mask. In this case, 
a path is an ordered sequence of chamfer mask vectors { V,^, Vjj, Vj,} with 
having its origin at p, 'V'i^ having its end at q, and ii,i2, ■■■,ii G [l...m]. 
In other words, Vj G [1../] Vj. G C and pc| = + + ... + Vj,. As Q 
is an Abelian group and a module on TZ, the order of the vectors does not 
matter, and we can bring together the different vectors ( of the path 
equal to the same mask vector V^, k G [l..m]. Observe that / is the number of 
vectors within the path and m is the size of the chamfer mask. Moreover, as 
for each vector G C, — G C, we can consider only positive coefficients. 
We obtain the following definition: 
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Definition 2.8 (Path from p to q) Given a chamfer mask 

C = {(/^ k,uik)ke[i..m] ^ Q X T^} and two points p, q G ^, a path Vp^ from p 

to q is a sequence of vectors of the mask C such that: 

m 

^pq = '^k^k = p4 With Wk e [l..m], ak G 7^+ 

k=l 

Definition 2.9 (Cost of a path) The cost W of a such a path Vpq is de- 
fined by: 

m 

>V(^pq) = T.»k-Wk 
k=l 

Since a mask C contains a basis of Q, and is symmetric, such a path always 
exists for any couple of points (p, q) with positive coefficients. 

Definition 2.10 (Weighted distance) A weighted distance dc associated 
with a chamfer mask C between two points p and ci in Q is the minimum 
of the costs associated with paths Vp^ linking p to q. 

(ic(p,q) =min{>V(Ppq)} 



2.3 Weighted distances properties 



Theorem 2.1 A weighted distance is invariant under translation. 

Proof. This proof can be found in [28]. Given any three points p, q, r in 
Q, dc{r + p,r + q) = dc{p, q). Indeed, (r + p)(r + q) = pc| and any path 

"^7 — , — \7 — , — t = ^k=i Oik^k = (r + p)(r + q) = p4 from r + p to r + q is 
(r + p)(r + q) 

also a path from p to q, with the same cost: W(V-, r-. 1) = VVfPf^) = 

(r + p)(r + q)' P^' 



Corollary 2.1 A property of a weighted distance between any two points (q, r) 
of Q can be expressed as a property of a weighted distance between the origin 
O and a point p. 

Taking p = r — q leads to dc{0, p) = dc{0 + q, r — q + q) = dc{q,, r). 

In the following, we will express the properties of weighted distances from the 
origin to any point p of ^ and denote (ic(p) = dciO, p). 

To be able to forecast the weighted distance map, from a single point, from 
the chamfer mask C , we divide C into several sectors spanning the points of Q 
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such that the distance to the origin of a point lying within a sector will only 
depend on a restricted number (actually n - the dimension) of mask weights. 
We define the following objects: 

Definition 2.11 (M-sector) Given a family of n independent vectors ofQ, 
(v^fc)fcg[i n]' ffi-sector ( Vi, V2, V„) is the region o/M" spanned by the 
vectors Vi, V2, i.e.: 

Definition 2.12 (^-sector) The ^-sector (( Vi, V2, is the set of 

points belonging to Q which are included in the M-sector (Vi, V2, v^„). 
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Fig. 1. Example of M-sector and t?-sector in (Z^,Z, +, •). 

Remark 2.2 It is possible that a ^-sector does not correspond to the set 

of points p such that p = J2k=i ^k^k, ^k G T^^- For example, as illustrated 

^ 

in Figure 1, if we take the module (Z , Z, +, ■), the point p(l, 1) lies in the M- 
sector (Vi = (2,1), V2 = (1,2)) ((1, ij = 1/3 ■ (2, ij + 1/3 • (1,2)) and as 
it is a point of Z^, it also lies in the ^-sector ((^1,^2)). But it cannot be 
written as Ai ■ Vi + A2 ■ v^2 with Ai, A2 G Z. ♦ 

Definition 2.13 (Wedge of a chamfer mask) We call a wedge of a cham- 
fer mask C, a Q -sector formed by a family T = (/^ijj^^i n ''^ vectors of C 
which does not contain any other vector of C. 

To avoid the situation illustrated in Remark 2.2 we consider only ^-sectors 
based on a basis of Q. 

Definition 2.14 (^-basis-sector) We call a Q-sector ((Vi, ... where 
i^k)k=i n ^'^ ^ basis of +, ■) a Q -basis-sector. 

By definition of a basis, a ^-basis-sector corresponds exactly to the the set of 
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points p such that 1? = Ylk=i -^fc^fc) -^fc G T^'^ ■ 

Given a family JF = {'^k)k£[i n] '^^ ^ independent vectors, we denote A^r G TZ 



A°. = det(Vi,V2,...,V„) 



vl vl ■■■vl 



Vi V2 ■ ■ ■ f 



and VA; G [l-.r?-], we consider the function : 



g — >n 

It I — > A^(l?) 



such that: 



A^(3?) = det(Vi, Vfc-i, 3?, Vfc+i, -^n) 



,,2 



■ ■ V 



Lemma 2.1 A family T = {'^k)ke[i n] ^'^ ^ &aszs 0/ {g,TZ, +, ■), zif 



e ^, VA; e [1.. 



Proof. The proof in (Z^, Z) can be found in [29]. 

By definition of a basis, a family JF = ( Vfc)^^^-^ is a basis of {Q,7l, +, ■) iff 
Vx G ^, 3ai, ^2, On G 7?. such that x = Y.k=i ■ '^k- 
As we consider sub-modules of M", x G M" and V/c G G M". If 

the vectors of JF are not independent, is not a basis of Q, and as A^r = 0, 
Wit G ^, ^A^(l?) ^ 7^. If JF is an independent family of vectors, as 

is a vector space, JF is a basis of (M",M) and 3ai, ^2, G M such that 

X = J2k=i oik ■ v^fc- This can be written: 



X 



yx J 



vl vl 



■ vt 



\ 



X 



\Vi ■■■ vlj 



«2 



M 
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As T is an independent family, A^r = det(Aljr) ^ and the matrix A^jr can 
be inverted such that 



' 0L\ 


\ 


' v\ 


v\. 




-1 




\ 


0L2 




v\ 


vl- 


■■vl 










/ 




■ 


■■vl) 






/ 



Inverting M.jr using Cramer's rule leads to 
VA; e [l..n], Ofc = X A^, 



and JF is a basis of {G,7V) iff VA; G [1..?^], at G 7?.. ■ 

When {Q,TZ) = (Z",Z) we obtain the condition A^r = ±1 for cone regularity 
defined in [30]. 

We can always organize C as a set of wedges (taking ra-tuples of vectors that do 
not contain any other vectors). To avoid the situation in Remark 2.2 and to be 
able to forecast the final weighted distance from the chamfer mask, we choose 
masks whose wedges are all ^-basis. When {Q, TZ) = (Z^, Z), taking every cou- 
ple (if, V) of adjacent vectors of C in clockwise (or counter-clockwise) order 
gives such an organization. If, for each couple det(if, V) = ±1, then, they 
are a ^-basis (cf. [24]). For example, in Figure 2 of Remark 2.3, {(( Vi, v^2>) , 
((^2,^3)), ((Va - Vl)), ((-Vl -Va}), (("^2-^3)), ((- ^3, Vl})} is an organi- 
zation of C in G-basis wedges. For n > 3, this organization may be more 
complicated. Indeed several ways of organizing n+1 independent vectors into 
two wedges may exist: for example, if we take the vectors Vi = (1,0,0), 
= (1,1,0), V3 = (1,1,1), V4 = (1,0,1) and their symmetric vec- 
tors, the wedges ((^1,^2,^3)) and {(^3,^4,^1)) are a ^-basis, but we 
can also consider the wedges (( v^4, v^i, ^^2)) and ((^2,^3,^4)). In [31], an 
automatic recursive method is given; it is based on Farey series to organize 
chamfer masks of (Z'^, Z) with ^-basis wedges. In the general case, considering 
a mask C containing only a ^-basis JF = ( Vfc, Wk)k=i..n and their symmetric 
wedges JFg = {—^k,Wk)k=i..n, leads to two ^-basis wedges (if Vl? E Q,Wk E 
[l..n],^A^(l?) G 7^, ^A%{lt) = ±^A^(1?) G 7^ as 7^ is a group). 

Moreover, the other wedges are JF; = ( Vi, — V^, V^+i, V„) for 
/ G [l..n]. They all are ^-basis wedges as V3? G Q,'il G [l..n], -^A^^ (x) = 

±-^A^(x) G 7?. as 7^ is a group. We can then add vectors to this mask, taking 
care of keeping an organization in ^-basis wedges. 

Remark 2.3 Given a chamfer mask C = {i^k,Wk)ke[i..n] E Q x 71} 
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which is organized in ^-basis wedges, and given any point p lying in a wedge 
(( Vj^, Vj^)) of C, there exist a path Vp = Z]fc=i «fc v^fc with Vfc G [l..n] 
(i.e. the point p can by reached by a linear combination of only n vectors 
among the m vectors of the mask C with coefficients in TZ~^). However, the 
final weighted distance at this point may not be a linear combination of the 
n weights corresponding to these n vectors in the chamfer mask. 



For example, in (Z^,Z), let us consider the mask C containing the weighted 
vectors (v*i,ioi) = ((1,0), 2), (v*2,-!i'2) = ((2,1), 5), (v*3,-!«3) = ((1, 1), 1) and their symmet- 
ric vectors. The families :Fi = (v*i, v*2) and :F2 = {^2,^3), generating the wedges 
5i = ((v*i,v*2)> and 52 = ((v*2,V3)) respectively, are a basis of (Z^,Z). Indeed, 

1 2 



AO 



1 



1 and VI? G Z2, 



2 G Z and: 



1 1 
X Al (It) = - X 



x^ 2 

X2 1 



1 ■ a;^ - 2 ■ a;^ G Z 



X 



1 x^ 




1 ■ - ■ G Z. 



In the same way for JF2, we have: A^r^ 
and: 



1 



xAUlt) 



X 



A^ 



1 

X2 1 

2 x^ 

1 X2 



2 1 
1 1 



1 and V3? G Z, x\ x^ G Z 



1 ■ x^ - 1 ■ x^ G Z 



2 ■ x^ - 1 ■ x^ G Z 



However, as illustrated in Figure 2 the point p = (2, 1) lying in the ^-basis- 
sector Si can be reached only by the path Vi = 1- Vi + 1 ■ v^2 = 1 ■ + 1 ■ v*i 
containing only vectors of J^i and VV(7'i) = wi+W2 = 2 + 5 = 7. But p can 
also be reached by the following path V2 = 2 ■ Vi + 1 ■ V3 = 1 ■ V3 + 2 • Vi 
belonging neither to the sector Si nor to the sector ^2 with a cost yV{V2) = 
2XW1 + W3 = 2x2 + 1 = 5. As the weighted distance of a point (with respect 
to the origin) is the minimum cost of all paths allowed by the mask, we have 
rfc(p) < 5 and thus rfc(p) ^ W(Pi). ♦ 
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J >WP(2,1) 



S. 



Fig. 2. Example where the weighted distance does not depend on ^-basis sector 
vectors. 



To avoid the situation mentioned in Remark 2.3, we add restrictions to the 
mask weights. These restrictions rely on the fact that the polytope formed by 
the chamfer mask vectors normaUzed by their weights is convex. 

Definition 2.15 (Normalized chamfer mask polytope) We call the poly- 
tope o/ M" whose faces are the n — 1- dimensional pyramids formed by the n 
vectors of each wedge of a chamfer mask Aic = {{'^k,Wk)ke[i..m\ G ^ x 71} 
normalized by their weights, i.e. Vj^. = — ■ Vj^, k G [l..n] for each wedge 

(( Vj^, Vj^)) of C the normalized chamfer mask polytope, denoted Be ■ 
Note that as {Q,TZ, +, ■) is a module but not a vector space, these points may 
not be in Q . 

Lemma 2.2 // the normalized polytope Be is convex, the weighted distance 
of any point p lying inside a wedge (( Vj^, ... Vj^)) of a chamfer mask C = 
{( Vfc, G ^ X 71} depends only on the weights Wi^, ...,Wi^. 

Proof. This proof can be found in [24] for (Z^, Z). 

If p = O, the proof is obvious. Let a point p 7^ O, p G (( V^^, Vj^)) be 
given. As JF = ( Vj^, Vj^) is a ^-basis, there exists a path P = J2k=i ^k'^i^ 
from O to p containing only vectors of having positive coefficients (V/c G 
[L.n], Ofc > 0). The cost of this path is W{V) = Y.k=i '^k x Wi^. We can write 
V as 



n 



dk X Wi^ 



— ■ V,, = WiV) ■ 



V = [^aiXWiA X 



E 



X 



1=1 



k=l 
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with 



u -p 



E 



Uk X Wi 



'Ik- 



Since ifp is a convex combination of the n normahzed vectors of JF: V/c G [1 •.■".], 
< v-^-"*'^"'''' — < 1 and y]t_i "''^"''fc — = 1 ligg Qn Moreover, as 

the faces of Be are convex {n — 1-dimensional polytope with n vertices) if-p 
also hes on the face formed by the family oi Eq-, i.e. on the boundary of Be- 
Let us now consider another path Q = X^fcLi Pk'^k from O to p containing 
arbitrary vectors of C. As C is symmetric, we can take V/c G [l..n],/5fc > 
without loss of generality. Then >V(Q) > W(P). Indeed, 

Q = (e ■ E >< ^ ■ = ms) ■ 

Vfc=l / k=l ^k=l Pk X Wk Wk 



with 



Q 



E 



Pk X Wfc 



Er=i Pk X 



1? Q is a convex combination of m normalized vectors of C, and lies thus within 
the convex polytope Bq. Moreover, we have Op = V = W{V) ■ l?-p = W(Q) ■ 
if Q = Q with l?p and 1? q having the same direction (Op) with positive 
coefficients. W-p \\'vtv\\ = Wq ||"i?q||. As if g lies within Be and 'vt-p lies on 
the boundary of Sc, ||lfp|| > ||"ifQ|| (as i3c is centered in O, if ||1?q|| > ||"ifp||, 
ug would be farther from the origin than a point of the boundary, and thus 
outside Be) thus W(P) < W(Q) and rfc(p) = W(P). ■ 
Figure 3 shows an example of vectors if-p and "ffg of a mask C G (Z^, Z). 

Corollary 2.2 // the vertices of each face of the convex hull conv(Sc) of Be 
are normalized vectors corresponding to Q -basis sectors of C, then the vectors 
of C whose corresponding normalized vectors do not lie on the convex hull of 
Be are not used to compute the final weighted distance. 

Proof The proof can be found in [32,28] for vectorial spaces. Suppose there 
exists a vector G C such that v ; = —V; does not lie on the con- 

vex hull of Be (i.e. lies within Be). Let us consider the ^-basis sector 
(( Vj^, Vj^)) in which lies and whose corresponding normalized vec- 
tors form a face of conv(i3c). There exists a path V = Ylk=i ^^k^ik ~ 
from O to V; and with VV(P) < wi. Indeed, V; = wi^i = W{V)~itp with 
if p = Y,k=i — V jj. a convex combination of n vectors of JF lying within 



13 



Fig. 3. Example where Be is convex. We consider the vectors 
vi(l,0), V2(l,l), V3(0,l), V4(-l,l) and the mask C = {(vi,4), (v2,5), (v3,4), 
(— vi, 4), (— V2, 5), (— V3, 4)}. We consider the path "P = vi + V2 and Q = 2 • vi + V3 
from O to p(2,l). We then have ||ug|| < ||up||. 

JF. We thus have v/ < and wi > W(V). Finally, as a weighted dis- 

tance is the minimum of the costs of all possible path, in any path containing 
V;, V/ will be replaced by V whose cost is lower. Note that this is what 
happens in Remark 2.3. ■ 



Remark 2-4 If there exist faces of conv(i3c) formed by normalized vectors 
whose corresponding mask vectors are not ^-basis, the vectors of C which does 
not lie on conv(i3c) may be used, but this leads to a final weighted distance 
which may not be homogeneous along some directions, and we also may not 
be able to forecast the weighted distance inside a wedge with the only vectors 
generating the wedge. 

For example, if we consider the mask C with the weighted vectors: Vi = 
((i;o5, 3), V2 = ((m], 2), V3 = ((Ojj, 3), V4 = ((-1,1), 2) and their sym- 
metric vectors, as shown in Figure 4, dc{0, (0, 2)) = 4 7^ 2 x dc{0, (0, 1)) = 
2x3. Moreover, dc{0, (0, 2)) does not only depend on V2 and V3 or only on 
V3 and V4. Note that the vectors ^2 and ^4 generating the corresponding 
face of conv(i3c) do not form a basis (det(Vi, V2) = 27^ il)- ^ 



In the following, we consider only chamfer masks whose normalized polytope 
is convex. Indeed, if this is not the case, the mask may be redundant (Corol- 
lary 2.2), or even if not, we may not be able to forecast the final weighted 
distance, (Remark 2.3). We can note that this condition implies that coUinear 
vectors (which are note opposite vectors) are removed from the mask. 
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f 1 \ 

-V2 -V3 ■V4 

Fig. 4. When conv(i3c) 7^ Be the weighted distance may not be homogeneous. 

Considering previous lemmas and remarks, we can re-define a chamfer mask 
with stronger conditions as follows: 

Definition 2.16 (Chamfer mask (restricted)) A Chamfer Mask Cp is a 
finite set of weighted vectors {{'^kiWk)ke[i..m] G ^ x TZ} which satisfies the 
following properties: 



{positive weights) \/k, Wk > and Vfc 7^ 0(1) 

{symmetry) {'^,w) G Cp =^ {—^,w) G Cp (2) 

Vp G ^ 1 " 

{Organized in Q -basis wedges) f P ~ '^k^ik (3) 

v',„)>eCp J k=i 

{Convex normalized polytope) = conv{Bcp) (4) 

Theorem 2.2 Given a chamfer mask Cp = ^{^k,'Wk)k£[i..n] G ^ x 7^|, de- 
fined as in Definition 2.16, the weighted distance of any point p lying in a 
wedge (( Vj^, ... Vj^)) can he expressed by: 

1 " 
k=i 



Proof. This formula was given for (Z^, Z) in [31] without the entire proof. Let 
p be a point of Q lying in the wedge W = (( v'j^, ...Vj^)) of Cp. As 1^ is a basis 
of {Q, TV) (condition 3 of Definition 2.16) there exists a path V from O to p and 
«!, ...,Q!fc G 7^^ such that V = Op = Z]fcLi«fc v^ife- The proof of Lemma 2.1 
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gives yk G [l..n],afc = x A^(p) and the proof of Lemma 2.2 gives that 
the weight of P is minimal Thus dcp(p) = W(P) = ^ x J2k=i ^^(p) x ^ik- ■ 



Theorem 2.3 ^4 weighted distance computed with a chamfer mask 

Cp = {( Vfc, Wfc)fcg[i m] ^ ^ X 71} as defined in Definition 2.16 is a norm on 

Proof. Here, we have to show that dc is definite, positive and symmetric and 
satisfies the triangular inequality and the positive homogeneity properties. By 
definition of a weighted distance (Definition 2.10), given any points p G ^, 
there exist ai, 0^2, a™ G such that dcp{0, p) = (icp(p) = YJk=i c^kWk and 
Pp = Op = p = ET=i oik^k (cf. Definitions 2.8, 2.9 and 2.10). 

(1) Positivity (Needs conditions 1 and 2 of Definition 2.16) 

By condition 2 of Definition 2.16 we can choose ak > for all k G [l..n]. 
By Condition 1 we also have Wk > for all k G [l..n]. We thus have 
dc,ip)>0. 

(2) Definitivity (Needs conditions 1 and 2 of Definition 2.16) 

If p = O, then the path V = J2T=o ' is a path from O to p and its 
cost is W(P) = Er=oO xwk = 0. Thus, rfcp(p) = min(W(p^)) < 0. By 
positivity, we also have (icp(p) > 0. It follows that o?Cp(p) = 0. 
Conversely, if (icp(p) = 0, as VA; G [l..m], Wk > 0, we have VA; G [l..m], 
ttfe = and ^ = Er=i ■ Vfc = O". Thus p = O. 

(3) Symmetry (Needs the condition 2 of Definition 2.16) 
By condition 2 of Definition 2.16, we have that 

m 

T'^ = T.M-^k) (6) 
^ fc=l 

is a path from p to O and ^('^^^) ~ ^'-'^Op''' ^^^^ '^'-'^^ 

is minimal. Indeed, let us consider a path Q— ^ = y]I!^i 0k^k such that 

pO 

W(Q^) < >V(P^). Then Q^j = ELi A-(-^fc) is a path from O to 
p and VV(Qq-^) < VV('Pq-^) = c?Cp(p) which is impossible by definition 
of a weighted distance (Definition 2.10). Thus (icp(p, O) = VV(P— ^) = 

(4) Triangular inequality (Needs conditions 1 and 2 of Definition 2.16) 
The proof can be found in [33] for (Z^, Z, +, ■) and in [34] (Corollary 3.4) 
for a general module. We want Vp, q G dcp{p,c\) < 'icp(p) + c^Cp(q)- 
Let "Pq = J2k=i Pk ■ v^fc be the minimum cost path between O and q 
(i.e. rfcp(q) = j:T=iPkWk), and by Eq. (6), P_p = Ek=iM-^k). As 
the mask is symmetric, V-p can be written V-p = YlT=i Ik'^k for some 
7fc and >V(P__p) = W('Pp). Let us suppose dcp{p,q) > rfcp(p) + c^Cp(g)- 
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The path Ppq = P_p + = Er=i(7fc + f3k)^k = p6 + = 
is a path from p to q. As (JZ, +) is an Abehan group, VV('Pp^) = 

W(P_p) + W(Pq) = rfcp(p) + c?Cp(q) < c^Cp (P,q) which is impossible 
by definition of a weighted distance (Definition 2.10). By contradiction, 
we have rfcp(p, q) < dc^{p) + dc^iq). 

Remark 2.5 Note that a decomposition in ^-basis wedges is not needed 
for this condition. The only condition needed is to be able to extract a 
basis among all mask vectors. This is the case, for example, for masks 
containing only vectors corresponding to knight displacements. Indeed, 
each wedge of this mask is not a Z-basis (see Figure 1 and Remark 2.2. 
However, this mask induces a metric, [35] . ♦ 

Remark 2.6 The triangular inequality does not depend on the choice 
of the weights. ♦ 

(5) Positive homogeneity (Needs conditions 1, 2, 3 and 4 of Definition 2.16) 
The proof can be found in [30] for (Z'^,Z, +, ■) and in [28] for a general 
module. Let be A G 7^. Let W = (( Vj^, Vj^)) be the wedges of Cp in 
which p lies. By Theorem 2.2, we have (icp(p) = x Z]fc=i ^^(p) ^ 'Wik- 
If A > 0, the point A ■ p also lies in the wedge W (it is a point of Q, and 
Ap = ELi A X cik^ik with yk e [l..n],X x a^, G 7^+). By Theorem 2.2 
we have: 

1 " 

dcM ■ P) = ^^(A ■ P) X ^u. 

k=l 

1 " 

k=i 
1 " 

= ^ E A X det (Vi,, ^, Vi^^^, ... X 

k=i 

1 " 

^.5^ fc=l 

= A X rfcp(p) 

If A < 0, (icp('^ ■ p) = dcp{—\ ■ p) with —A > (symmetry property) 
and thus, dc^X ■ p) = -A x dcp{p) = \M x dc^ip)- ■ 



2.4 Weight optimization 



In addition to their metric and norm properties, weighted distances can be 
made more invariant to rotation. As a weighted distance is obtained by com- 
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puting the smallest weight of several paths between two points, the first im- 
provement to obtain a weighted distance with high rotational invariance is 
to allow a larger number of allowed directions for the paths. This means in- 
creasing precision by increasing the number of weighted vectors of the chamfer 
mask. 

To increase accuracy, another way is to choose suitable weights for mask vec- 
tors. This more challenging issue as been often addressed in the literature. 
The first optimal chamfer weights computation was performed for a 2-D 3x3 
mask in a square grid [25]. Then authors computed optimal weights with dif- 
ferent optimality criteria [36,31], for larger masks [37,36,38] and for anisotropic 
grid [39,40,41]. Authors of [41,31,42] proposed an automatic computation of 
optimal chamfer weights for rectangular grids. 

Observe that computing the distance transform and finding optimal weights 
is not directly related to the problem of estimating the length of straight lines 
in a discrete image [43]. For optimal weights for such estimations, see [44]. 

In all of the previous papers, the computation of optimal chamfer weights is 
performed the same way: 

(1) First, a chamfer mask is built and decomposed in wedges. 

(2) Then, the final weighted distance from the origin to an arbitrary point of 
the grid is expressed. The variables corresponding to the mask weights are 
unknown, but variables corresponding to vector coordinates are known. 
In the general case, given a chamfer mask defined as Definition 2.16, and 
a point p E Q lying in the wedge JF = (( Vj^, Vjj, ... Vj^)) the value of 
the weighted distance at this point is 



(3) In the same way, the Euclidean distance from the origin to this point 
p is expressed. In the general case, given a grid which is a sub-module 
of M"" with an elongation of Si, S2, s„ in each canonical direction, the 
Euclidean distance between the origin and the point p can be expressed 
in the following way: 



(4) After these three steps, the error between the weighed distance and the 
Euclidean one can be expressed for any point p G ^. This error can be 
either absolute (the difference between these two values) [37] or relative 
(the difference is divided by the Euclidean distance) [36]. The error can 
be expressed as follows: 



1 



■J^ k=i 
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E{p)=dc^{p) -dEip) 



\k=i 




n 



The general relative error can be expressed as follows: 



Erelip) 



dcpjp) - dEjp) 
dsip) 



ELl A^(P) X 



- 1 



(5) The previous errors are n-dimensional functions of the coordinates p^jp"^, 
... p"- of p. To reduce the number of these variables and be able to find 
extrema, the maximal error is computed either on a hyperplane or on a 
sphere. For example, in three dimensions, the error can be computed on 
a plane X = T,Y = T or Z = T [37,31] or on the sphere of radius T [42]. 
This error function is continuous on a compact set (the n — 1 pyramid 
formed by the n vectors of each wedge). Thus it is bounded and attains 
its bounds. These bounds can be located either at the vertices of the 
pyramid or inside (including other bounds such as edges) the pyramid. 

(6) Computing the derivatives of the error function gives the point Pmax at 
which the error is maximum (due to the sign of the derivatives) inside the 
wedge. In this way, the maximum error within the wedge Emax = -E(pmax) 
can be obtained (if Pmax lies within the wedge). 

(7) The other extrema (minima) Ei, E2, En are obtained for the n vectors 
delimiting the wedge (the vertices of the n—1 pyramid). When computing 
the error on a sphere of radius R, in the general case, the extrema can 
be expressed in the following way (V/ G [l..n]): 



with II V^ll = yJ2k=ii^kViy being the Euclidean norm of the vector V/ 

expressed in the world coordinate. 

(8) Minimizing the maximum error leads to computing optimal real weights 
with the following equation: Emax = —Ei = —E2 = ... = —En. 

(9) A depth-first search in an integer weights tree taking the error into ac- 
count can lead to optimal integer weights set for a given mask [31,42]. 

We have generalized weighted distance properties found in the literature to 
modules. These properties are true for the well-known cubic grid, but can 
also be applied to other grids such as the FCC and BCC grids as we will 
see in Section 4. In applications, an efficient algorithm for computing distance 
transforms is is needed. In the following section, we discuss such an algorithm: 
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the two-scan chamfer algorithm. We prove that the algorithm produces correct 
result for images on general point lattices. 



3 How to compute distance maps using the chamfer algorithm 

There are basically three families of algorithms for computing weighted dis- 
tance transforms - bucket-sort (also known as wave-front propagation), paral- 
lel, and sequential algorithms. Initially in the bucket-sort algorithm [45], the 
border points of the object are stored in a list. These points are updated with 
the distance to the background. The distances are propagated by removing 
the updated points from the list and adding the neighbors of these points to 
the list. This is iterated until the list is empty. The parallel algorithm [46] is 
the most intuitive one. Given the original image, a new image is obtained by 
applying the chamfer mask simultaneously to each point of the image giving 
the minimum distance value at each of these points. This process is applied 
to the new image. By applying the procedure iteratively until stability, the 
distance map is obtained (a proof using infimal convolution is found in [34]). 

Rosenfeld & Pfaltz [12] showed that this is equivalent to two sequential scans 
of the image for grids in 1? and a 3 x 3 mask. Sequential means that we use 
previously updated values of the same image to obtain new updated values. 
Two advantages are that a two-scan algorithm is enough to construct a dis- 
tance map in any dimension and that the complexity is known (0(M), where 
M is the number of points in the image). The worst case complexity of the 
parallel algorithm is O(M^) - M iterations of operations on M points. This 
result is now generalized and proved to be correct on any grid for which all 
grid point coordinates can be expressed by a basis. 

In this section, 5 is a finite subset of Q and 7?. = Z. The cardinality of S 
is denoted M = card(5). Observe that the chamfer weights are in TZ~^ = N 
(Definition 2.16, condition 1) and that the weighted distance have values in 
7^+ (Definition 2.10). 

Definition 3.1 Given a = (ai, 02, ... , an), where oi, 02, . . . , ctn ^ and a G 
M. Let 

'^a,u = {p e ^ : ai'p^ + 02/ + . . . + anP^Ua] , 
where □ is one of the relations <, <, =, >, or >. 

For example, T^^ is a hyperplane and T^<, '^a>i ^'^d 7^% are half-spaces 
separated by the hyperplane T^^- 

Definition 3.2 (Scanning masks) Let a chamfer maskC be given. Let also 
ai, a2, ■ ■ ■ yttn E TZ defining the hyperplane Tq_ such that TQ_r\C = (/) be given. 
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The scanning masks with respect to C are defined as 

In the following, the notation Vj G Ci and G C^, / G {1, 2}, means that the 
pair ( G Ci will be used. 

Definition 3.3 (Scanning order) A scanning order is an ordering of the 
M = card(iS) points in S, denoted pi, p2, . . . , Pm- 

For a scanning mask to propagate distances efficiently, it is important that, 
in each step of the propagation, the values at the points in S from which the 
mask can propagate distances have already been updated. This is guaranteed 
if each point that can be reached by the scanning mask either has already 
been visited or is outside the image. 

Definition 3.4 (Mask supporting a scanning order) Let pi, p2, . . . , Pm 

be a scanning order and Ci a scanning mask. The scanning mask Ci supports 
the scanning order if 

Vpi, G Ci, 3i' <i : (pj/ =Pi + or pi + ^ S) . 

Proposition 3.1 Given a scanning maskCi and an image S, there is a scan- 
ning order such that Ci supports the scanning order. 

Proof. Let 7jL such that T^^^ fl C = be the hyperplane defining Ci and C2 
(we consider Ci here, the proof for C2 is similar). Now we consider two sets, 
iS\ V and SUV, where S is the complement of S and the elements in V are the 
already visited points in S. Let V = {pi, p2, . . . , Pn}, N < card(5). For every 
p E S, there is a dp G M (iTp = aip^ + a2P^ + . . . + a„p") such that p G X^^-- 
Since 5 \ V is finite, there is a point p G 5 \ V with corresponding ap such 
that p = {S\V)nTf^^^, S\V = {S\V)nT^l^, and = {S\V)nT^l^. Thus, 
W E Ci, ai ■ p^ + a2 ■ p"^ + ■ ■ ■ + an ■ p"" = cTp and ai-v^ + a2-v'^ + . . . + an-v"' < 0, 
so ai ■ {p^ + v^) + a2 ■ (p^ + u^) + . . . + an- {p"' + v"-) < Up and thus p + V G ^p,< 
which has an empty intersection with iS \ V, so p + V ^ {S \V). Let this p 
be Pat+i, assign it to V and repeat until S = V. M 

Definition 3.5 (Chamfer algorithm) Let the scanning masks Ci and C2 
and scanning orders pi, p2, . . . , Pm and pm, Pm-i, ■ ■ ■ ,Pi such that the masks 
support these scanning orders be given. 

• Initially, Vp G X, /(p) < hoo and Vp G X, /(p) ^ 0. 

• The image is scanned two times using the two scanning orders. 
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For each visited point pi in scan I, 
I 

f{pi) ^ mill 



mm 



J 



Remark 3.1 Usually, in the first scan f{pi) is omitted, as f{Pi) is only 
computed for X, so we know /(j) = oo. ♦ 



Now we prove that, in a distance map DMx, for any point p' on a path 
defining the distance between p G X and its closest background point q G X, 
the value in the distance map at point p' is given by DMx{p') = d{ci, p'). In 
other words, the closest point in the background from p' is q. 

Lemma 3.1 Let p G X and q G X &e such that d{p,q) = d{p,X) and 
I = {i : {'^ti,Wi) G C}. There is a set {ai G TZ} such that the point p can be 
written p = q + ^ v^j and d{p, q) = ^ ctiWi. 

For any set {Pi E TZ : < Pi < ai}, ifp' = q + Yj f^i^i ^ ^^^^^ d{p',X) = 

c^(p',q) = Y^PiWi- 

iei 

Proof. Let d(j)',X^ = K ^ J2 Pi^i- Since J2 Pi^i is the length of a path 
between p' = q + 5] Pi'^i and q, the weighted distance can not be larger than 
this. Assume that K < J^PiWi- Then there is a q' G X such that 

K = d{c^,p') <d (q, p') and 

d{'i',p)>d{q,p) = rf(p,x) . 

Now, 

d{q, P) = Y1 ^i^i = + ~ A)^* > d (q', p') + "^iai - Pi)wi. 



Since J2{<^i — Pi)wi is the length of one path (we cannot assume that it is a 
shortest path) between the points p' = q + Pi'^i and p = q + c^i'^i, it 
follows that Z(«i — Pi)wi > d (p', p) and thus 

d (q', p') + - > d (q, p) + d (p , p) > ci(q', p). 



which contradicts d (q', p) > d (q, p). ■ 
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To assure that the propagation does not depend on points outside the image, 
the definition of border points or of wedge preserving images defined below 
can be used. 

Definition 3.6 (Border point) Given a chamfer mask C, a point p E S is 
a border point if 

3V G C : p + V ^ 5. 

Lemma 3.2 If for all border points r G iS, r G X, then for a// p G X, q G X 
such that d{p,X) = d{p,q), all points in any shortest path between p and q 
are in S. 

Proof. Assume that a point p' in a shortest path between p and q is not in S. 
Then, since all border points are in the background, there is a background grid 
point in all paths between p and p'. Since d{p,p') < d{p,q), it follows that 
there must be a border point p" G X such that d{p, p") < d{p, p') < d{p, q) 
which contradicts the assumption d{p,q) = d{p,X). ■ 

Definition 3.7 (Wedge-preserving half-space) Given a chamfer mask C 
and a G M, ai, 02, • • • , ctn ^ such that 3p G 5 : p G X^=- The half-spaces 
T^^ and T^^ are wedge-preserving if 

Vp G r^%yW, either (W G W^, p + V G 7;'^<) or (W G W^, p + V G T^^^^) , 
where W denotes wedges in C. 

Definition 3.8 (Wedge-preserving image) Given a chamfer mask C, the 
image S is wedge-preserving if it is the intersection of wedge-preserving half- 
spaces. 

Lemma 3.3 // an image S is wedge-preserving, then for all p, q G S, all 
points in any shortest path between p and q are in S. 

Proof. Let S be wedge-preserving and let p, q G iS such that p = q + ^ «j 

with 2 = {i : G W}, where is a wedge. Assume that there is a fixed 
k eX and a set {7^ : < 7^ < and 7^ < a^} such that p' = q+^ li^i G S 

but p' + ^ S. In other words, if there is a point not in S in the shortest path 
between p and q, then some point p' in S has a neighbor p' -|- that is not in 
S. Since S is wedge-preserving, there is a plane defined by cr, ai, a2, . . . , a„ G 7^ 
such that p' G and p' + G T^^^. Since 

aip'^ + a2p'^ -|- . . . -|- a„p'" < a and 
aiij)'^ + vl) + a2(p'^ + vl) + . . . + ttnip'"' + ffc) > cr, we get 
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Olffc + 02^fc + . . . + Ctnffc > 0. 

Taking any r G 7^*^=, r + G ^'^>, which by Definition 3.8 impHes that 
Vr G Tj^=, W G VT, r + V G Tj^>. This imphes that 

y^eW,^e%^^. (7) 

Now, p = p' + ^iai- G 5, but 

= ai (p'' + ^1) + «2 (p" + V^,) + . . . + a„ (p'" + ^l) + 
+ ^ (ai-7i) (aiv] + asi;,^ + . . . + a„<) + (afc-7fc-l) (aivl + agf^ + . . . + a^i^^) > 

> cr + ^ (ttj - 7i)0 + (Ofc - 7fe - 1)0 = cr, 

since p' + G 7j^>, ctj > 7j if i 7^ A; and > 7fc + 1, and using Eq. 7. We 
have p G i.e. p ^ 5. Contradiction. ■ 

Remark 3.2 As the image generally is given first, a mask can often be 
built with preserved wedges and then cut into two scanning masks. Then the 
points in the image can be ordered such that the masks support the scanning 
directions and then the distance map can be computed using our algorithm. ♦ 

Theorem 3.1 If either 

• for all border points p G iS, p G X or 

• S is wedge-preserving, 

then the chamfer algorithm in Definition 3.5 produces distance maps as defined 
in Definition 2.6. 

Proof. Let p G X and q G X be such that (i(p, q) = (i(p, X) and X = 
{i : (v^j,u;j) G C}. There is a set {a^ G TZ} such that the point p can be 
written p = q + ^ Cii^i and (i(p, q) = X! 

Let / be the image after the first scan. By Lemma 3.1, using that the masks 
support the scanning order, the local distances are propagated in the following 
way. Also, by Lemma 3.2 and Lemma 3.3, all points below are in S. 

/(q + VfcJ = Wk-^ for any ki such that 

Vfci G Ci, and > 1 
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sequential masks 

(a) 































J 




1 








00 


00 


00 





00 


3 


2 


1 





1 


00 


00 


00 


00 


00 


5 


4 


3 


2 


1 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 


00 















(b) 



(c) 



00 


00 


00 


00 


00 










00 


00 


00 





00 








00 


00 


00 


OO 


00 








00 


00 


00 


00 


00 










00 


00 


00 


00 


00 



4 


3 


2 


1 


2 










3 


2 


1 





1 








4 


3 


2 


1 


2 








5 


4 


3 


2 


3 










6 


5 


4 


3 


4 



(d) 
























































00 


00 











1 


1 











00 


00 


00 








1 


2 


1 








00 


00 


00 








1 


1 


1 
















































(f) 



(g) 



Fig. 5. Consider the mask C = {((-1, 0), 1), ((-1, 1), 1), ((1, 0), 1), ((1, -1), 1)} 
(a). The image in (b) is not wedge-preserving when using this mask. In (c), 
the result of the chamfer algorithm using Ci = {((—1,0), !),((— 1,1),!)} and 
C2 = {((1,0), 1), ((1, —1), 1)} is shown. By considering the wedge-preserving image 
in (d), the correct distance map (shown in (e)) is produced. By assigning all border 
points to the background, the image in (f) is achieved and the chamfer algorithm 
produces the correct result shown in (g). 



/(q+ v'fc, + V 



fcj) = Wfci + for any ki, k2 such that 

and Ofc^ , > 1 
(afci > 2 if fci = /C2), 



/|q+ 5Z "fc,^fc. 
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Let g be the image after the second scan. Using the notation q^v = q + 

(^k,^k,, we get 

for any li such that 

e C2, and ai^ > 1 
for any li, I2 such 
that V/^, G C2 
and a^^, > 1 
{ai, > 2 if /i = /s) 



This proves that all steps of a shortest path between p and q have been 
processed. Indeed, after the second scan, the value of the point p is (^(p) = 

c^(p,q) = Y^oiiWi. ■ 

Remark 3.3 In Figure 3, the need for the conditions in Theorem 3.1 are 
shown in example images. ♦ 



4 Best weights for the BCC and FCC grids 

The BCC grid (B) and the FCC grid (F) are defined as follows: 
B = {p(x, y, z) G and x = y = z (mod 2)} and 
F = {p(a;, y, z) e I? and x + y + z = {) (mod 2)}. 

A voxel is defined as the Voronoi region of a point p in a grid. Figure 6 (a) 
shows a voxel of a BCC grid. A BCC voxel has two kinds of face-neighbors 
(but no edge- or vertex-neighbors), which results in the 8-neighborhood (Fig- 
ure 6 (b)) and the 14- neighborhood (Figure 6 (c)). On the FCC grid, each voxel 
(see Figure 6 (d)) has 12 face- neighbors and 6 vertex-neighbors. The result- 
ing 12- and 18-neighborhoods are shown in (Figure 6 (e)) and (Figure 6 (f)), 
respectively. 

Note that the high number of face neighbors (12 for the FCC grid and 14 
for the BCC grid) implies these grids are more compact than the cubic grid. 



fi'(qjv + V^/i) = 9{f\N)+Wi^ 

g{(iN + v^/i + v^/2) = g{<iN) + wi^ + wi^ 



qN+ (^k^h\ =9i(lN) + Y ^h^k 
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(a) (b) (c) (d) (e) (f) 

Fig. 6. Different neighborhoods for the BCC grid ((a)-(c)) and the FCC grid 
((d)-(f)). 

which has only 6 face neighbors. 

4.1 Results of the previous sections applied to BCC grids 

To apply the results of Sections 2.3 and 3 to the BCC grid, we should check 
that: 

• B is a sub-module of 

• Masks decomposed in B-basis wedges can be created 

• A scanning order can be defined on the BCC lattice. 

Lemma 4.1 LetM"' = {p : p^,p'^, ...,p"- E X and = p"^ = ... = p"^ (mod 2)} 
, then (B'",Z, +, x) is a sub-module o/M". 

Proof. Let us prove first that B" is an Abelian group. Given p, q G B"-, 
there exist 02, •••,an, and P2, ■■■,Pn ^ ^ such that Wk e [2..n] p^ = p^ + 2ak 
and = + 2/3^. The point r = p + q = q + p also belongs to B". Indeed, 
Wk G [2..n],r'' =p^ + q^ = p^ + 2ak + q^ + 2pk = ri + 27fc with 7^ = afc+/9fc e Z. 
Moreover, — p G B" as —p^ = —p^ + 2{—ak) with —Ok G Z. We now have to 
prove that VA G Z, Vp G B'^, A ■ p G B*^. Indeed, VA; G [l..n], \p^ G Z and 
VA; G [2..n]\p^ = Xp^ + 2(Aafc) with (Aofc) G Z. ■ 
This result is true for all n G N \ {0} and obviously for n = 3. 

4.1.1 Chamfer mask geometry 

Lemma 4.2 A family T = ( Vi, v^2, v^a) is a basis o/B iff 
A°. = det(Vi,V2,V3) = ±4. 

Proof. First, we prove that any determinant of 3 vectors Vi, V2, v^3 of B is 
a multiple of 4. Indeed, given a vector '^ti{xi,yi, Zi), G B iff there exists 
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ctj, /3i G Z such that i/i = Xi + 2ai and Zi = Xi + 2(3i. Then the determinant of 
Vi, and V3 is 



det(Vi, ^2, V: 





Xi X2 X3 






X2 


2^3 


= 


yi 2/2 2/3 




Xi + lax 


X2 + 2^2 


X3 + 2^3 




^1 ^2 ^3 




xi + 2/5i 


X2 + 2/^2 


X3 + 2/33 




Xi X2 X3 








= 4 


ai 0L1 03 










/?1 /?2 /53 










e z 









With the notations of paragraph 2.3 we obtain for any vector y, z) 
{x,x + 2a,x + 2p) of B: 



AO 



1 



Xi X2 X3 
«! 0:2 «3 
A /52 /93 



X 4 



X X2 X3 

a 0^2 ^3 

/3/32 /?3 



and 4r x A^(l?) G Z iff 



±1 which means A^r = ±4. The same 



Xi X2 X3 
ttl 0^2 Q^3 
/?1 /52 /53 

result apphes for x A^(lf) and x A^(lf). By Lemma 2.1 we obtain 
is a basis of B iff A^r = ±4. ■ 



To compute a weighted distance map on an image stored on a BCC grid, one 
can consider a mask built using 8-neighbors: 



C8 = {(V,(±l,±l,±l),ti;i),} 



This mask contains only one type of weight as the Euclidean distance between 
each face-sharing neighbor to the central voxel is the same. 
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This mask can be decomposed into B-basis sectors. Indeed, let us consider 
the wedge jFg = 1, 1), V2(— 1, 1, 1), ^^3(1, —1, 1))), and its symmetric 

wedges. 



A° =det(Vi,V2,V; 



1 -1 1 
1 1 -1 
1 1 1 



and from Lemma 4.2, is a B-basis sector. 

If we want larger masks, we can split the wedge and its symmetric wedges 
according to the Farey triangulation technique [31,30]. To do so, given a B- 
basis wedge JF = ((^1,^2,^3)), choose an edge to be split, say the edge 
between Vi and V2. We create V(a;, z) such that V = Vi © ^2 defined 
by X = xi + X2, y = yi + 2/2 and z = zi + Z2. The two obtained wedges 
= (( ^2, v^s)) and = (( v^i, ^5 ^3)) ^^^e B-basis sectors. Indeed, 



AO 



det(V, V2, 



= ±4. 



Xi + X2 X2 X3 

vi + y2 y2 vs 

Z\ + Z2 Z2 Zz 





Xi X2 X3 




X2 X2 X3 




yi 2/2 2/3 


+ 


2/2 2/2 2/3 




Zl Z2 Z3 




Z2 Z2 Z3 



In the same way, Ajr^ = ±4. 

If we apply the first steps of this method, we can retrieve the chamfer mask 
using the 14-neighborhood proposed by [27]. But we can also build larger 
masks by splitting the new wedges. 



4.1.2 Chamfer mask weights 

To compute optimal integer chamfer weights for BCC grid the depth-first 
search method of [31,42]^ can be used. This method computes weights sets 
leading to weighted distance maps with small relative error with respect to 
the corresponding Euclidean map. Moreover, it keeps a weights set only if 
the normalized polytope of the corresponding mask is convex, which allows 
application of the results from Section 2.3. 

The following tables give sets of integer chamfer mask weights for BCC grids 
with the scale factor allowing comparison of the final weighted distance map 

^ The corresponding code is available at www.cb.uu.se/~tcl8/code-data-set 



29 



with an Euclidean one composed with real numbers. They also give the max- 
imum relative error that can occur between the weighted distance map and 
the Euclidean one. 

For the two first cases, we also give the optimal real weights computed using 
the equations of Section 2.4. 

Balls obtained by using different number of weights (shown in bold below) are 
shown in Figure 7. 

(1) One weight 

For the mask corresponding to the 8-neighborhood, using only one weight, 
the equations of section 2.4 give the following result: 



real optimal weight 


1.268 


real optimal error (%) 


26.79 


integer optimal weight 


1 


real map scale factor 


1.268 



These results are consistent with [27]. 
(2) Two weights 



(3) 



vector 


weight 


real 








weig 


hts 










(1 1 1) 


Wl 


1.547 


1 


2 


3 


4 


5 


6 


13 


19 




(2 0) 


W2 


1.786 


2 


3 


4 


5 


6 


7 


15 


22 




Scale factor 


1 


1.268 


0.731 


0.504 


0.383 


0.308 


0.256 


0.119 


0.081 




Error (%) 


10.69 


26.79 


15.59 


12.70 


11.60 


11.07 


10.78 


10.72 


10.71 




Three 


weights 




















vector 


weight 










weight! 










(1 1 1) 


Wl 


1 


2 


4 


5 


6 


13 


19 


26 


33 


(2 0) 


W2 


2 


2 


5 


6 


7 


15 


22 


30 


38 


(2 2 0) 


W3 


2 


3 


7 


8 


10 


22 


31 


43 


54 


Scale factor 


1.268 


0.899 


0.396 


0.325 


0.270 


0.125 


0.0857 


0.0626 


0.0494 


Error (%) 


26.79 


10.10 


8.50 


7.94 


6.39 


6.34 


6.12 


6.12 


6.11 










Optimal 


error 


(%): 6.02 











Note: for the previous masks, we displayed real optimal weights to be 
able to compare our results with previous papers. However, as real op- 
timal weights are computed sector by sector, they may change from one 
sector to another for the same mask vector. We do not display them for 
several sectors as they may not be consistent. 



(4) Four weights 
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vector 


weight 


weights 


(1 1 1) 


Wl 


1 


2 


4 


5 


6 


9 


15 


26 


(2 0) 


W2 


2 


2 


4 


6 


7 


10 


17 


29 


(2 2 0) 


W3 


2 


3 


6 


8 


10 


14 


24 


41 


(3 1 1) 


W4 


3 


4 


7 


10 


12 


17 


29 


50 


Scale factor 


1.268 


0.899 


0.460 


0.334 


0.275 


0.194 


0.113 


0.0662 


Error (%) 


26.79 


10.10 


7.94 


5.57 


4.73 


4.21 


4.00 


3.99 








Optimal error 


{%): 


3.96 









4-2 Results of previous sections applied to FCC grids 

To apply the results of Sections 2.3 and 3 to the FCC grid, we should again 
check that: 

• F is a sub-module of 

• Masks decomposed in F-basis wedges can be created 

• A scanning order can be defined on the FCC lattice. 

The calculations are similar to those in Section 4.1. 

Lemma 4.3 Let F" = {p : p^.p^, G Z and p^ + p^ + p^ = Q (mod 2)} 
, then (F",Z, +,-) is a sub-module o/M". 

Lemma 4.4 A family JF = (Vi, V3) is a basis of¥ iff 

A0, = det(Vi,V2,V3) = ±2. 

These Lemmas can be proved similar to Lemma 4.1 and 4.2 by noting that 
for any p G F", there exist a G Z such that p^ = p^ + . . . + p""^ + 2a. 

To compute the weighted distance on FCC, we first consider the smallest mask, 
i.e., the mask containing 12- neighbors: 

C12 = {(V,(±l, ±1, 0), Wl), (V,(±l, 0, ±1), w,), (V,(0, ±1, ±1), Wl), } 

The wedges = ((Vi(l, 1,0), ^2(1, 0, 1), ^3(1, -1, 0))), 

J-'f = (( ^1(1, 1, 0), ^2(1, 0, 1), ^4(2, 0, 0))), and their symmetric wedges are 

F-basis sectors. 

= = 2 

■^F F 

and from Lemma 4.4, jFp and JF| are F-basis sectors. 
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The splitting of the sectors are analogous to Section 4.1.1, but the determi- 
nants equals ±2 instead of ±4. 



We get the following mask weights. Balls obtained by using different number 
of weights (shown in bold below) are shown in Figure 7. 



(1) One weight 

For the mask corresponding to the 8-neighborhood, using only one weight, 
the equations of section 2.4 give the following result: 



real optimal weight 


1.172 


real optimal error (%) 


17.16 


integer optimal weight 


1 


real map scale factor 


1.172 



These results are consistent with [27]. 
(2) Two weights 



(3) 



vector 


weight 


real 


weights 












(1 1 0) 


Wl 


1.271 
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1 
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(2 0) 


W2 


1.798 
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Scale factor 
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1.464 


1.172 


0.636 












Error (%) 


10.10 


26.79 


17.16 


10.10 












Three 


weights 
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(2 0) 


W2 


1 


2 


3 
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10 


16 


22 


(2 1 1) 
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7 


10 


12 


19 


26 


Scale factor 


1.464 


1.172 


0.694 


0.636 


0.325 


0.226 


0.191 


0.121 


0.0887 


Error (%) 


26.79 


17.16 


15.04 


10.10 


7.94 


7.76 


6.19 


6.16 


5.95 










Optimal 


error (%): 5.93 











(4) Four weights 
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ghts 








(1 1 0) 


Wl 
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3 
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12 


(2 0) 
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7 


7 


13 


17 


(2 1 1) 


W3 
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2 
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5 


9 


9 


16 


21 


(2 2 2) 


W4 


2 


3 


5 


7 


12 


13 


23 


30 


Scale factor 


1.268 


1.172 


0.651 


0.472 


0.274 


0.272 


0.150 


0.113 


Error (%) 


26.79 


17.16 


7.94 


5.57 


5.15 


4.64 


4.63 


4.07 








Optimal error 


(%): 3.98 
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5 Conclusions 



We have presented a general theory for weighted distances. This allows appli- 
cation of the weighted distance transform to any point-lattice. The optimal 
weight calculation and the construction of a chamfer algorithm to produce 
distance maps are straight-forward by the procedure in Section 2.4 and The- 
orem 3.1. For the chamfer algorithm to construct correct distance-maps, the 
only limitation is that the image must be wedge-preserving, Definition 3.7, or 
have all border points in the background. Definition 3.6. With these condi- 
tions which can easily be satisfied, images in any dimension can be considered. 
It is also worth mentioning that despite its popularity, to our knowledge, no 
proof of the correctness of the two-scan algorithm has been published until 
now, except for the 2D 3x3 mask algorithm proof from 1966 [12]. Also, the 
important general formula for the weighted distance between two grid points 
in Theorem 2.2 is new. 

The application of these properties and of the chamfer algorithm are numerous 
for the cubic and parallelepiped grids in the literature. To illustrate how this 
theory can be applied. Section 4 gives examples for the FCC and BCC grids. 
But it can be applied to any point-lattice, for example the 3D grid where the 
voxels are hexagonal cylinders as suggested in [47] or 4D a grid as in [48], or 
even other grids as long as they satisfy the module conditions. 
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